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ABSTRACT 


A  computer  program  for  calculation  of  the  echo  area  of 
smoothly  joined,  N  section  convex  conducting  surfaces  of  revolution^ 
described  by  a  second  degree  equation  is  presented.  For  the  case  of  E  q 
(parallel)  polarization  of  the  incident  and  scattered  fields  the 
solution  is  obtained  by  a  combination  of  geometrical  optics  and 
creeping  wave  theory.  For  the  case  of  E^  (perpendicular)  polari¬ 
zation  the  solution  is  obtained  using  geometrical  optics,  and  the 
creeping  wave  is  neglected.  The  computed  results  for  Eq  polari¬ 
zation  are  in  good  agreement  with  measurements  for  prolate 
spheroids,  prolate  spheroid-sphere,  and  prolate  spheroid-oblate 
spheroid  combinations. 
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A  COMPUTER  PROGRAM  FOR  BACKSCATTER 
BY  SMOOTHLY  JOINED,  SECOND  DEGREE 
SURFACES  OF  REVOLUTION 


I.  INTRODUCTION 

The  computer  program  listed  in  Appendix  1,  applies  the  geo  - 
metrical  optics  and  creeping  wave  solutions  described  in  Ref.  1,  2 
to  obtain  the  backscattered  fields  of  a  surface  of  revolution.  The 
target  is  composed  of  N  sections,  each  section  described  by  a 
second  degree  equation.  In  addition  the  target  must  be  convex,  and 
be  smoothly  joined  at  the  boundaries  between  sections.  The  program 
as  listed  computes  the  backscattered  field  and  echo  area  for  Eq 
(parallel)  polarization  as  a  function  of  the  incidence  angle,  for  a 
given  wavelength.  The  geometrical  optics  backscattered  field  is 
independent  of  polarization,  thus  the  geometrical  optics  scattered 
field  for  the  case  of  E^  (perpendicular)  polarization  is  also  obtained. 
The  creeping  wave  scattered  field  has  not  been  included  for  E  ^ 
polarization  due  to  the  difficulty  in  obtaining  the  ray  path  of  the 
creeping  wave  in  this  case.  1  This  program  has  been  tested  for 
E  g  polarization  for  prolate  spheroid,  prolate  spheroid-sphere, 
and  prolate  spheroid-oblate  spheroid  targets.  The  results  of  these 
tests  are  presented  in  Ref.  1. 

II.  TARGET  DESCRIPTION 

The  surface  is  described  in  each  section  by  the  second 
degree  equation 

(1)  F(r,6 )  =  Ax  r2  sin2  0  +  A2  r2  cos20+A3r2  sin0cosO 

+  A4  r  cos  9  +  A5  r  sin  0  +  A6  =  0 

The  constants  .  .  .  A^  specify  the  surface  in  the  section  (i)  bound¬ 
ed  by  the  angles  0^  and  0^  +  i .  The  surface  which  may  be  represent¬ 
ed  by  Eq.  (1)  includes  the  sphere,  prolate  and  oblate  spheroids,  and 
the  ogive.  In  general  any  surface  derived  from  a  conic  section  can 
be  described  by  this  equation*  More  complex  surfaces  may  be 
represented  by  using  a  large  number  of  sections  and  approximating 
the  desired  surface  by  a  second  degree  surface  within  each  section. 
This  program  provides  for  20  sections  but  this  number  can  be  readily 
increased. 
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A  restriction  on  the  target  specification,  which  is  a  result  of 
the  method  for  finding  the  specular  point,  is  that  the  coordinate  origin 
be  located  such  that  the  normal  to  the  surface  at  z  =  0  be  r  directed. 
That  is,  dr/d0|z  =  q  =  0  .  Discontinuities  in  the  derivatives  (dnr/dGP) 
at  the  junctions  between  sections  may  exist.  However,  the  effects  of 
such  discontinuities  as  well  as  the  effects  of  tips  are  not  included  in 
this  program.  A  discontinuity  in  the  first  derivative  (dr/dQ  at  the 
junction  (nwedge"  discontinuities)  causes  diffraction,  and  can  be 
evaluated  using  wedge  diffraction  techniques.  l>  2  A  discontinuity 
in  the  second  derivative  (d2  r/d^P  )  at  the  junction  results  in  reflect¬ 
ion  and  transmission  of  an  incident  creeping  wave  at  such  a  discon¬ 
tinuity.  An  example  of  this  effect,  the  spherically  capped  ogive,  has 
been  discussed  in  Refs.  1  and  2.  As  the  effects  of  a  discontinuity  in 
the  first  derivative  are  significent  this  program  may  not  give  good 
results  for  such  a  target.  The  creeping  wave  reflection  effect  due 
to  a  moderate  second  derivative  discontinuity  may  be  neglected  with 
good  results.  1  Thus  good  results  for  the  scattered  field  of  perfectly 
conducting  convex,  smoothly  joined  surfaces  of  revolution  may  be 
obtained  using  this  program. 

III.  THE  CREEPING  WAVE  COMPUTER  PROGRAM 

The  creeping  wave  computer  program  given  in  Appendix  I  uses 
geometrical  optics  and  creeping  wave  theory  to  calculate  the  back- 
scattered  field  of  the  target.  The  geometrical  optics  field  is  obtain¬ 
ed  by  identifying  the  specular  point  and  calculating  the  Gaussian 
curvature  at  the  specular  point.  The  backscattered  field  is  then 
calculated  as  given  in  Ref.  1.  The  creeping  wave  backscattered 
field  is  obtained  by  identifying  the  points  of  attachment  and  re¬ 
radiation  of  the  creeping  wave,  calculating  the  diffraction  coefficient 
at  these  points,  and  by  computing  the  attenuation  of  the  creeping  wave 
on  a  path  defined  by  the  E-plane  of  the  target.  The  backscattered 
field  is  then  calculated  as  given  in  Ref.  1.  These  contributions  are 
added  to  obtain  the  total  backscattered  field. 

Referring  to  the  computer  program  listing  shown  in  Appendix  I  ? 
the  function  of  the  significant  sections  of  the  program  will  be  dis¬ 
cussed.  The  card  numbers  associated  with  each  section  will  be 
specified.  This  discussion,  together  with  the  comment  cards  in¬ 
cluded  in  the  program  listing,  is  intended  to  give  sufficient  infor¬ 
mation  about  the  program  to  enable  a  qualified  programmer  to  both 
use  and  modify  the  program.  Statements  which  are  in  common  use 
in  Fortran  IV,  such  as  DIMENSION  ,  COMPLEX,  and  FORMAT 
statements  will  not  be  discussed  as  it  is  assumed  that  the  reader 
has  a  knowledge  of  Fortran  IV. 
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The  COMMON  declaration  (0011)  is  used  to  store  the  constants 
required  in  Eq.  (1)  in  the  common  block  labelled  DATA/.  This 
common  block  is  used  in  conjunction  with  the  unlabelled  common  block 
to  transfer  a  particular  set  of  constants  A1  (I)  to  A6  (I)  into  the  un¬ 
labelled  common  regions  shared  by  the  subroutines.  This  provision 
reduces  the  number  of  calling  variables  required  by  each  subroutine. 

The  COMMON  declaration  (0013)  is  used  to  store  and  manipu¬ 
late  the  angles  corresponding  to  the  section  boundaries. 

The  statements  0018-0022  initialize  constants  which  are  re¬ 
quired  in  the  calculations. 

The  READ  statements  0023-0027  read  the  required  data, 
and  provision  is  also  made  to  write  out  this  data  for  the  purpose  of 
identification. 

The  statements  0028-0030  set  up  the  incrementation  of  the 
incidence  angle  THT  and  the  statement  00  31  calculates  the  propaga¬ 
tion  factor  FK. 

Next  the  loop  for  incrementing  the  incidence  angle  is  entered, 
and  the  angle  converted  to  radians.  The  logic  statements  (0036,  0037) 
serve  two  purposes  if  the  incidence  angle  exceeds  90°.  First  the 
incidence  angle  is  constrained  to  the  range  0.  <  THT  <  tt  2.  Then 
the  subroutine  FNVRT  (N)  is  called.  The  subroutine  FNVRT  (N) 
inverts  the  target,  i.  e.  ,  the  portion  of  the  target  defined  as  section 
#1  becomes  section  #N ,  section  #2  becomes  section  #(N-1)  and  so 
forth.  This  provision  is  necessary  so  that  the  propagation  direction 
of  the  creeping  waves  which  start  at  the  point  of  attachment  is  always 
the  positive  angular  direction.  This  provision  simplifies  the  logic 
required  to  compute  the  creeping  wave  contributions .  Figure  1 
illustrates  this  provision. 

The  propagation  vector  and  polarization  vector  of  the  inci¬ 
dent  wave  are  computed  next  (0038-0046).  These  vectors  are  used 
in  the  determination  of  the  specular  point  and  the  points  of  attach¬ 
ment  and  reradiation  of  the  creeping  wave. 

Next,  the  location  of  the  specular  point  is  determined 
(0047-0060).  This  is  done  by  incrementing  by  DTHB  along  the 
surface  of  the  target  in  the  cj>  =  0  plane,  performing  the  scalar 
product  of  the  propagation  vector  and  the  surface  normal  vector  at 
each  point,  and  by  finding  the  point  at  which  this  product  is  a  maximum. 
This  point,  in  section  NSP,  is  the  specular  point  (THSPP). 
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Fig.  1.  A  target  (a)  and  the  inversion  (b) . 
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The  location  of  the  points  of  attachment  and  reradiation  of  the 
creeping  wave  for  $  =  0,  tt  is  determined  next  (0062-0090).  This  is 
done  in  the  same  manner  as  the  determination  of  the  specular  point 
location  except  that  the  scalar  product  of  the  polarization  and  normal 
vectors  is  used.  After  these  points  have  been  found  they  are  written 
out  together  with  the  specular  point  (0091,  0092). 

Having  determined  the  location  of  the  specular  point  (NSP, 
THSPP,  RSP)  the  Gaussian  curvature  at  the  specular  point  is  comput¬ 
ed  and  used  to  calculate  the  geometrical  optics  scattered  field  (0093- 
0103). 


The  creeping  wave  path  length  is  computed  in  0104-0129. 

This  is  done  by  starting  at  the  attachment  point  (THCWL,  <j>  =  0)  and 
performing  a  numerical  integration  of  the  differential  arc  length  along 
the  surface  to  the  point  6  =  tt  .  The  length  thus  determined  is  CWLl, 
This  process  is  then  repeated  for  the  attachment  point  (THCWL,  c|)  =  Tr) 
with  a  result  CWL2.  The  product  of  the  propagation  factor  and  the 
sum  of  the  path  lengths  then  gives  the  free  space  phase  of  the  creep¬ 
ing  wave  which  propagates  around  the  target. 

Next  the  complex  attenuation  of  the  creeping  wave  is  com¬ 
puted  (0130-0150).  The  same  process  used  to  compute  the  path 
length  is  applied,  with  the  exception  that  the  complex  attenuation 
coefficient  must  be  integrated.  This  is  accomplished  using  the 
EXTERNAL,  ALPHDS  as  the  integration  function. 

Having  obtained  the  creeping  wave  path  length,  attenuation, 
and  points  of  attachment  and  reradiation,  it  is  easy  to  calculate  the 
creeping  wave  scattered  field.  This  is  done  in  0151-0161.  The  phase 
is  first  calculated  and  the  square  of  the  diffraction  coefficient  de¬ 
termined.  Then  the  creeping  wave  scattered  field  ECW  is  computed. 

The  total  scattered  field  and  echo  area  in  square  wavelengths 
are  computed  in  0162-0165.  Next  the  common  regions  are  reset  if 
THTD  >  90°  and  the  incidence  angle  in  degrees,  the  geometrical 
optics,  creeping  wave,  total  scattered  field  and  echo  area  are 
written.  At  the  completion  of  the  loop  (0169)  the  program  is  ter¬ 
minated. 
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IV.  FUNCTIONS  AND  SUBROUTINES 


In  addition  to  the  arguments  required  to  call  each  subroutine 
or  function  as  given  in  the  description  below,  it  is  necessary  to  trans 
fer  the  constants  in  Eq.  (1)  for  the  section  (I)  from  the  /DATA  / 
common  block  to  the  unlabelled  common  block.  All  function  and 
subroutines  in  which  these  COMMON  statements  appear  require  that 
such  a  transfer  be  made.  This  provision  reduces  the  number  of  argu 
ments  required  in  the  functions  and  subroutines. 

Complex  Function  DELSP(THT)  0172-0179 

This  function  calculates  the  incremental  arc  length  in  the  0 
direction  at  THT.  Although  the  arc  length  is  a  real  number  this 
function  is  declared  complex  so  that  it  may  be  used  as  an  external 
function  in  the  numerical  integration. 

Subroutine  FNVRT(N)  0180-0195 

This  subroutine  interchanges  the  data  which  specifies  the 
target,  thus  inverting  the  target  by  180  in  0. 

Complex  Function  SPHASE  (THTI,  PHII  ,  THTB,  PHIB,  RB,  FK) 
0196-0204 

This  function  determines  the  phase  of  the  incident  field 
(THTI,  PHII)  at  a  point  on  the  target  (RB,  THTB,  PHIB)  for  a 
propagation  factor  FK. 

Complex  Function  PHASE  (THTI,  PHII,  THTB,  PHIB,  RB ,  FK) 

0205-0214 

This  function  determines  the  backscattered  phase  of  the  field 
at  (RB,  THTB,  PHIB)  for  an  incident  field  (THTI,  PHII)  and  prop¬ 
agation  factor  FK. 

Complex  Function  ALPH  (RH01,  RH02,  WAVE)  0215-0230 

This  function  computes  the  complex  creeping  wave  attenu¬ 
ation  coefficient  for  the  orthogonal  radii  of  curvature  RH01  and  RH02 
and  a  wavelength  WAVE.  RH01  is  the  radius  of  curvature  in  the 
propagation  direction. 

Complex  Function  DSQ  (RI1,  RSI,  WAVE)  0231-0243 

This  function  computes  the  square  of  the  creeping  wave 
diffraction  coefficient  as  a  function  of  the  radii  of  curvature  RI1, 

RSI  in  the  propagation  direction  and  the  wavelength  WAVE. 
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Complex  Function  ALPHDS  (THT)  0244-0260 

This  function  determines  the  product  of  the  creeping  wave 
attenuation  coefficient  and  the  metric  of  the  surface  which  is  needed 
in  the  integration  for  the  attenuation  along  the  path. 

Subroutine  DIFFGO  (R ,  THT,  RTHT,  RTTH ,  ECAP,  FCAP,  GCAP, 

ELC  ,  F LC  ,  GLC)  0261-0301 

This  subroutine  calculates  the  first  (RTHT)  and  second  (RTTH) 
derivatives  of  the  distance  (R)  from  the  origin  to  the  surface  at  the 
point  (R,  THT).  In  addition  the  coefficients  of  the  first  and  second 
fundamental  forms  of  differential  geometry  of  the  surface  ECAP, 
FCAP,  GCAP,  EEC,  FLC,  GLC  are  computed.1 

Subroutine  FSPDT  (RSP,  THSP,  FNT,  THTS ,  THTF ,  DTHT  ,  PHI, 

VX,  VY ,  VZ)  0302-0326 

This  subroutine  finds  the  point  (RSP,  THSP)  in  the  section 
FNT  bounded  by  THTS  and  THTF  where  the  scalar  product  of  the 
surface  normal  vector  and  the  vector  (VX,  VY,  VZ)  is  a  maximum. 
This  is  done  by  incrementing  in  angle  by  DTHT  and  selecting  the 
largest  scalar  product. 

Subroutine  FINT  (SSS,  FCTI,  FLL,  FUL,  ERRR ,  NX)  0327-0377 

This  is  a  subroutine  for  numerical  integration  of  the  complex 
external  function  FCTI  between  the  limits  FLL  and  FUL.  The  com¬ 
plex  result  is  returned  in  SSS.  ERRR  specifies  the  percent  error 
desired  and  the  integer  NX  determines  whether  equal  integration 
increments  (NX  =  1)  or  adjusted  increments  (NX  =  2)  are  used.  A 
description  of  this  integration  technique  is  given  in  Ref.  3. 

Subroutine  FNORM  (FNVX,  FN VY  ,  FNVZ,  R,  THT,  PHI)  0378-0400 
The  subroutine  FNORM  calculates  the  surface  normal  vector 
(FNVX,  FNVY,  FNVZ)  at  the  point  on  the  surface  described  by  the 
spherical  coordinates  R,  THT,  PHI. 

Function  RAD(THT)  0401-0419 

This  function  computes  the  distance  from  the  origin 
to  the  point  on  the  surface  at  the  angle  THT. 

Subroutine  FCOMM  (I)  0420-0431 

This  subroutine  transfers  the  constants  required  by  Eq.  ( 1) 
for  the  I^1  section  from  the  /DATA/  common  storage  block  to  the 
unlabelled  common  block. 


7 


V.  DATA 


A  set  of  input  data  is  shown  in  Fig.  2,  The  order  of  cards  is 
as  follows : 

Card  1  -  specifies  the  number  of  sections  (N)  in  the  target. 
Card  2  -  specifies  the  wavelength  and  the  increment  in 
incidence  angle. 

Next  N  cards  -  specify  the  initial  and  final  angular  boundaries 
of  the  section  in  radians  and  the  constants 
required  in  Eq.  (1). 

The  particular  target  specified  by  this  data  is  the  prolate 
spheroid-oblate  spheroid  combination  described  in  Ref.  1. 


00002 

1,0  2.0 

0.  1.570796  2.4799  .5102  0.  0. 

1.570796  3.1415927  2.4799  5.0955  0.  0. 


Fig.  2  .  The  input  data. 

VI.  CONC  FUSION 

A  computer  program  for  backscatter  by  smoothly  joined, 
second  degree  surfaces  of  revolution  has  been  developed  using  the 
theory  presented  in  Refs.  1  and  2.  This  program  has  been  tested 
for  prolate  spheroids,  prolate  spheroid-sphere,  and  prolate  spheroid- 
oblate  spheroid  combinations  with  good  results.  1  The  description  of 
the  computer  program  given  in  this  report  should  enable  a  capable 
programmer  to  use  and  modify  this  program. 

The  creeping  wave  computer  program  was  originally  in¬ 
tended  to  be  integrated  with  the  computer  program  based  upon  wedge 
diffraction  techniques.  1  ’  4  This  goal  can  be  obtained  by  modification 
of  the  logic  sections  of  the  two  programs,  and  by  converting  the 
computational  sections  of  each  of  the  programs  to  a  subroutine  form. 
Thus  a  program  utilizing  both  wedge  diffraction  and  creeping  wave 
theory  can  be  assembled. 
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APPENDIX  I 

THE  COMPUTER  PROGRAM  LISTING 
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SEXECUTE  I B JOB  0000 

$ I B JOB  GO. MAP  0001 

SIBFTC  CREEP  LIST. NODECK  0002 

1  FORMAT (115)  0003 

2  FORMAT (2F10. 5 )  0004 

3  FORMAT (8F 10.5 )  0005 

4  FORMAT (3F15.8)  0006 

5  FORMAT (7F15. 8 )  0007 

6  F ORM AT ( 7H  THSPP  =  F 1 5 . 8 . GH  THCWU  =  F 1 5 . 8 . 7H  THCWL  =  F15.8)  0008 

7  F0RMAT(5H  NSP= 115. 6H  NSWU=115.6H  NCWL=115)  0009 

8  F0RMAT(6H  CWL  1=F15»8, 6H  CWL2  =  F15.8)  0010 

COMMONRA 1  .RA3.RB1  . RA9 . RA 1 0 . R A  1  1  . WAVE/DA TA/AR 1  ( 20 ) . AR3 (20 ) «BR1 (20 ) ♦  0011 

CAR9 ( 20 ) , AR1 0 ( 20 ) , ARl 1 (20 )  0012 

C0MM0N/ANGL/THT1 (20) .ThTF(20),ThTP(20)  0013 

COMPLEX  PALPI. PALP2  0014 

COMPLEX  PHASP.EGO. ALPHDS.ATTEN.ALPl .ALP2.PHCW1 . PHCW2 .DSOC.ECW. ETOT  0015 
COMPLEX  CPTHL .PHASE . ALPH .DSQ . DELSP.CML 1 .CML2. CEXP  0016 

COMPLEX  SPHASE  0017 

DEGRAD=0. 01745329  0018 

RADEG=57. 29578  0019 

PI=3. 1415927  0020 

PI2=PI/2.  0021 

TP=2.*PI  0022 

READ ( 5 . 1 )  N  0023 

READ (5.2 )WAVE  «DTHT  0024 

WRI  TE  (6. 2  MAVE.DTHT  0025 

Read (5.3 ) (THT I < IR) ♦ THTF (  IR ) , ARl  <  IR  )  , AR3  <  IR ) . BRl  ( I R  > . AR9(  IR) . AR10 ( I  0026 

CR) .ARl  1  ( IR)  .  IR=l  .N)  0027 

DT  HB  =  DEGRAD  002Q 

FNT= 1 80./DTHT  0029 

NT  =  FNT- 1  •  0030 

FK=  TP/WAVE  0031 

DO  100  NTH  -  1  . NT . 1  0032 

FNTH=NTH  0033 

THTD=FNTH*DTHT  0034 

THT=DEGRAD*THTD  0035 

IF ( THTD.GT.90.  )  THT  =P 1  —  THT  0036 

IF (THTD.GT.90. )  CALL  FNVRT(N)  0037 

C  CALCULATE  INCIDENCE  VECTOR  0038 

PHI =0 •  0039 

VIX=-SIN(THT)  0040 

V IZ=— COS(THT )  0041 

VIV=0.  0042 

C  CALCULATE  INCIDENT  E-VECTOR-PARALLEL  POL.  0043 

E I X  =  COS ( THT  )  0044 

E I  V  =  0 .  0045 

E I Z=-S I N ( THT )  0046 

C  DETERMINE  SPECULAR  POINT  0047 

FNSPP=0.  0048 

NSP=0  0049 

RSPP  =  0 .  0050 

THSPP=0 •  0051 

D04 1 0  ISP=I.N«1  0052 

CALL  FCOMM(ISP)  0053 

CALL  F SPOT (RSP. THSP.FNSP. THT 1 ( I SP ) .THTF (ISP) . DTHB. 0. .VlX.VIY.VIZ)  0054 

I F (FNSPP.GT .FNSP.OR. THSP.GT .P I  2  )GOTO  411  0055 

FNSPP= ABS (FnSP  )  0056 

RSPP=RSP  0057 

THSPP=THSP  0058 

NSP=ISP  0059 

411  CONTINUE  0060 

410  CONTINUE  0061 

C  DETERMINE  CREEPING  WAVE  POINTS  0062 

FCWU=0.  0063 
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51  I 

510 

52  1 

520 

C 

2  1 

20 

22 

C 

522 

524 


0064 

0065 

0066 

0067 

0060 

0069 

0070 

0071 

0072 

0073 

0074 

0075 

0076 

0077 

0078 

0079 

0080 

0001 

0002 

0003 

0004 

0005 

0006 

0087 

0000 

0009 

0090 

0091 

0092 

0093 

0094 

0095 

0096 

0097 

0090 

0099 

0100 

0101 

0102 

0103 

0104 

0105 

0106 

0  107 

0108 

0109 

0110 

0111 

0112 

0113 

0114 

0  115 

0116 

0117 

0110 

0119 

0120 

0121 

0122 

0123 

0124 

0125 

0126 

0127 


NCWU=0 

RSCWU=0. 

THCWU=0 . 

DO  51 0  I CW=I  iNtl 
CALL  FCOMM ( I CW ) 

CALL  F  SPOT ( RCWU  »  THWU » FWU »  THT !  ( I C W ) ,  THTF (  ICW ) »  DTHB  ,0«.EIX,EIY,EIZ  ) 

IF (FCWU.GT. FWU.OR.THWU.LT. PI  2 )G0  TO  51  I 

FCWU=ABS(FWU) 

RSCWU=RCWU 
THCWU=THWU 
NCWU  = ICW 
CONTINUE 
CONT I NUE 
FCWL=0. 

NCWL=0 
R  SCWL  =  O • 

THCWL  =  0 • 

DO  520  ICW=I iNi 1 
CALL  FCOMM (ICW) 

CALL  FSPDT (RCWL • THWL.FWL.THT,  ICICW).THTF(ICW). DTHB. PI  .EIX.EIY.EIZ) 

I F (FCWL.GT .FWL.OR. THWL.GT.PI 2 )  GO  TO  521 

FCWL=FWL 

RSCWL  =RCWL 

THCWL=THWL 

NCWL= I CW 

CONTINUE 

CONT I NUE 

WR I TE ( 6  «  7 )  NSP.NCWU.NCWL 
WRI TE (6.6 ) THSPP . THCWU. THCWL 
CALCULATE  GEOMETRICAL  OPTICS  TERM 
CALL  FCOMM ( NSP ) 

CALL  DIFFGO (RSPP. THSPP. RDUM.RDDM, EC AP.F CAP, GCAP.ELC. FLC .GLC ) 

IF(GCAP)  20,21.20 

GAUSS=ELC/ECAP 

GAUSS= GAUSS* GAUSS 

GO  TO  22 

G  AUSS= (ELC*GLC-FLC*FLC ) / ( EC AP*GC AP-FC AP*FC AP ) 

CONT I NUE 

PH ASP= PHASE ( THT .0. . THSPP . 0 . . RSPP.FK ) 

EGO=-SQRT ( 1 . /GAUSS ) *PHASP/2 . 

CALCULATE  CREEPING  WAVE  PATH  LENGTH 
EXTERNAL  DELSP 
CWL1 =0. 

DO  522  NCPU=NCWU.N, I 
CALL  FCOMM(NCPU) 

TCWI =THCWU 
T  CW2  =  T  HTF ( NCPU ) 

I F ( NCPU.GT.NCWU )  TCW1=THTI (NCPU) 

CALL  F I  NT ( CML 1 .DELSP .TCWI  ,TCW2,5.0,2) 

RCML1=REAL ( CML 1  ) 

CWL1=CWLI+RCML1 

CONTINUE 

CWL2  =  0 . 

DO  524  NCPL=NCWL ,N . 1 
CALL  FCOMM (NCPL) 

TCW 1 = THCWL 
TCW2  =  THTF (NCPU ) 

IF (NCPL.GT.NCWL )  TCW 1 = THT I ( NCPL ) 

CALL  F I  NT (CML2, DELSP. TCWI  . TCW2.5.0.2 ) 

RCML2  =RE AL ( CML2  ) 

CWL2=CLW2+RCML2 

CONTINUE 

WRITE (6. 8)  CWL1.CWL2 
FKL 1 =FK*CWL 1 
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FKL2=FK*CWL2  01 

FKLCW=FKLI +FKL2  01 

C  calculate  CREEPING  WAVE  ATTENUATION  01 

EXTERNAL  alphds  01 

ALPI = ( 0 • *  0 •  )  01 

DO  530  NCPU=NCWU*N* I  01 

CALL  FCOMM(NCPU)  01 

TC W I  - T HCWU  01 

TCW2=THTF (NCPU)  01 

IF(NCPU*GT.NCWU)  TCW1=THTI (NCPU)  01 

CALL  FINT(PALPl tALPHDSfTCWl *TCW2*5.0*2 )  01 

ALP  I =ALPI +PALP1  01 

530  CONTINUE  01 

ALP2= (0. *0. )  01 

00  532  NCPL=NCWL*N* I  01 

CALL  FCOMM(NCPL)  01 

TCWI=THCWL  01 

TCW2  =  THTF (NCPL  )  01 

IF(NCPL.GT*NCWL)  TCWI=THTI (NCPL)  01 

CALL  FINT(PALP2<ALPHDS*TCW1 *TCW2, 5*0*2 )  01 

ALP2 - ALP2+P ALP2  '  01 

532  CONTINUE  01 

ATTEN=ALP1 +ALP2  01 

C  CALCULATE  CREEPING  WAVE  01 

PHCWI =SPHASE ( TH T  *  0 •  • THC  WU  t  0 •  ♦ RSCWU  t  FK )  01 

PHCW2=SPHASE ( THT .0. ♦THOfc'L^O. ♦ RSCWL ♦ FK )  01 

CALLDIFFGO(RSCWU%THCWU»RDUM»RDDM,ECAPl *FCAP1 • GCAP 1 tELCI *FLC1 *GLCI )  01 

CALLDIFFGO (RSCWL f  THCWLfRDUM* ROOM, EC AP2tFCAP2t GCAP2  *ELC2  *FLC2  *GLC2  >  0  1 

FKAP1 =ELCI /ECAPI  '  01 

FKAP2=ELC2/ECAP2  01 

RHCWI =1 ./FKAPI  01 

RHCW2  = 1  • /FK AP2  01 

DSQC=DSQ (RHCW 1 *RHCW2*WAVE)  01 

E C W  =  —2 • *DSOC*PHCW I *PHCW2*CEXP ( -ATTEN+ (0* ♦  -  I  •  )*F<LCW  )  0  I 

ETOT=EGO+ECW  01 

EMAG=CABS (ETOT )  01 

S  I  GM A  =  2 • *TP*EMAG*EMAG  01 

S  IGMAL= 1 0.*ALOG 1 0 ( SIGMA )  01 

I  F ( THTD.GT .90 •  )  CALL  FNVRT(N)  01 

WR I TE ( 6*5 ) THTD. EGO t ECW t ETOT  01 

WRITE (6. 4 )THTD* SIGMA *SIGMAL  01 

l60  CONTINUE  01 

STOP  01 

END  01 

SIBFTC  DELSP.  LIST  01 

COMPLEX  FUNCTION  DELSP ( THT )  01 

R  =  R AD ( THT )  01 

CALL  D IFFGO ( R  «  THT  % RDUMt  RDDM ♦ EC AP ♦ FCAP ♦ GC AP *ELC , FLC ♦ GLC )  01 

ECAP=ABS(ECAP )  01 

DELSP=CMPLX(SQRT(ECAP) *0. )  01 

RETURN  01 

END  01 

ilBFTC  FNVRT.  LIST  01 

SUBROUTINE  FNVRT(N)  01 

COMMONRA 1  ♦  R  A3  «  RB I  * R A 9 ♦ R A  1 0  t R A  1  1  ♦ WA VE/DA TA/AR 1  (20 )  ♦ AR3 (20 ) ♦ BR 1  <20 ) ♦  01 

C AR9 (20)  .ARI0(20) ♦ ARl  I  (20)  01 

C  OMMON/ ANGL/TH  T I  ( 20 )  * THTF ( 20 ) • THTP( 20 )  01 

P I =  3 • 1 4 1 5927  01 

DO  I  0  I  =  1  tN* I  01 

BRI  (  I  ) =  — BR 1  (  I  )  01 

AR9 (  I  )  =  — AR9 (  I  )  01 

THTP (  I  )  =P I  — THTF (  I  )  01 

THTF ( I ) =P I —THT I ( I )  01 

THT 1(1)  =  TH TP (  I  )  01 


28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 


12 


10  CONTINUE  0192 

RETURN  0193 

END  0 1 94 

SIBFTC  SPHAS  LIST  0195 

COMPLEX  FUNCTION  SPHASE ( THT I  .  PH  I  I  ,  THTB , PH  I B , RB , FK )  0196 

FS= 1 .  0197 

TEST=ABS(PHI I-PHIB)  0198 

IF(TEST.EQ.O. )  FS=-I.  0199 

FL=RB*COS(THTI+FS*THTB)  0200 

FKL=FK*FL  0201 

SPHASE =CMPLX ( COS (FKL ) .SIN (FKL  )  )  0202 

RETURN  0203 

END  0204 

$  I BF  T C  PHAS.  LIST  0205 

COMPLEX  FUNCTION  PHASE ( THT I  »  PH  I  I  1THT81 PH  I B , RB  »FK 1  0206 

FS= 1 •  0207 

TEST=ABS(PHI I-PHIB )  0208 

IF ( TEST.EQ.O. )  FS=-I.  0209 

FL  =  2.*RB*C0S< THT1+FS*THTB  )  0210 

FKL  =FK*FL  0211 

PHASE  =  CMPLX ( COS ( FKL  )  , SIN (FKL  )  )  0212 

RETURN  0213 

END  0214 

SIBFTC  ALPH.  LIST  0215 

COMPLEX  FUNCTION  ALPH ( RHO 1 . RH02 . WAVE )  0216 

COMPLEX  ALPH, EX  0217 

EX=CMPLX(0. 86603, 0.5)  0218 

I F ( RH02 .EQ.C. )  GO  TO  10  0219 

RA=RHO 1 /RHO 2  0220 

U2= ( 1 .48/EXP ( 0.84+RA ) 1+0.20  0221 

GO  TO  1 1  0222 

10  U2=0 . 20  0223 

11  AR 1 =2 • * ALOG ( RHO 1 1/3.  0224 

FR1=1 ./EXP< ARl )  0225 

WV  1 = ALOG (WAVE  1/3.  0226 

FWV= 1 ./EXP < WV 1 1  0227 

ALPH=U2*FR1 *FWV*EX  0228 

RETURN  0229 

END  0230 

SIBFTC  DSQ .  LIST  0231 

COMPLEX  FUNCTION  DSQ ( R I  1  . RS 1  , WAVE  1  0232 

COMPLEX  DSQ, EX  0233 

EX  =  CMPLX (  .96593.-. 25882  1  0234 

U1 =0.270  0235 

A=SQRT ( R I 1*RS1 1  0236 

AL  =  ALOG ( A ) /3.  0237 

WV  =  2. *ALOG < WAVE  1/3.  0238 

F  A  =  EXP ( AL  )  0239 

FW=EXP(WV)  0240 

DSQ=U1*FA*FW*EX  0241 

RETURN  0242 

END  0243 

SIBFTC  ALPDS.  LIST  0244 

COMPLEX  FUNCTION  ALPHDS ( THT )  0245 

COMMON  ARl .AR3.BR1 ,AR9,AR10, ARl 1 .WAVE  0246 

COMPLEX  ALPHDS, ALP, ALPH  0247 

PI=3. 1415927  0248 

R  =R AD (THT)  0249 

CALL  DI FFGO ( R.THT, RTHT .RDDM. ECAP ,FCAP. GCAP.ELC ,FLC ,GLC 1  0250 

FMET=SQRT (ECAP)  0251 

RHO 1 =  EC AP/ELC  0252 

RH02=GC AP/GLC  0253 

IF ( THT.EQ.O. .OR.THT.EQ.PI )  RH02=RH0l  0254 

RHO 1 =ABS ( RHO 1 )  0255 
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RH02=ABS (RH02 )  0256 

ALP=ALPH (RHO 1 ,RH02 » WAVE )  0257 

ALPHDS=ALP*FMET  0258 

RETURN  0259 

END  0260 

$ 1BFTC  D1FG0  LIST  0261 

SUBROUTINE  D 1 FFGO (R*THT  «  RTHT  .RTTH.ECAP.  FCAP  < GCAP »  ELC  tFLCiGLC  )  0262 

COMMON  A1 . A3.B 1 . A9 . A  1 0 . A1  1 . WAVE  0263 

C  DIFFERENTIAL  GEOMETRY  PROPERTIES  OF  QUADRIC  SURFACE  OF  REVOLUTION  0264 

C  F (R.THT ) =A1 *R*R*SI N (THT ) *S 1 N (THT ) +A3*R*R*C0S ( THT ) *COS ( TH T ) +B 1 *R*R*  0265 

C  S IN (THT )*COS ( THT )+A9*R*C0S (THT ) +A 1 0*R*S IN ( THT )+Al 1 =0  0266 

C  THETA-PHI  COORDINATES  0267 

C  UPPER  AND  LOWER  CASE  F  =  0  0268 

C  INPUT  A1 .A3.B1 .A9.A10.A1 1 .R. THT  0269 

PHI =  0 •  0270 

ST  =  S I N ( THT  )  0271 

CT  =  COS ( THT  )  0272 

S2T  =  S IN (2.*THT  )  0273 

C2T=C0S <2.*THT )  0274 

ST2  =  ST  *ST  0275 

CT2=CT*CT  0276 

SP  =  S I N ( PH  I  )  0277 

CP  =  COS ( PH  I  )  0278 

U=A1*ST2+A3*CT2+B1*ST*CT  0279 

V=A9*CT+A10*ST  0280 

UTHT  =  A  1 *S2T-A3*S2T  +  B 1  * ( CT2-ST2 )  0281 

VTHT=-A9*ST+A 1 0*CT  0282 

UTTH=2.* ( A1 -A3 )*C2T-2.*B1 *S2T  0283 

VTTH=— A9*CT-A 1 0*ST  '  0284 

Sl=2.*R*U+V  0285 

S2=R*UTHT+VTHT  0286 

FMAG=SQRT (SI *Sl +S2*S2 )  0287 

FMAG2=FMAG*FMAG  0288 

FMAG3=FMAG2*FMAG  0289 

C  CALCULATE  UPPER  CASE  E » F « G  0290 

RTHT  =  — R*S2/S 1  0291 

ECAP  =  R*R+RTHT*RTHT  0292 

FCAP  *  0.  0293 

GCAP  =  R*R*ST2  0294 

C  CALCULATE  LOWER  CASE  E.F.G  0295 

RTTH=-(2 .*PTHT* (2. *R*UTHT+VTHT )+2.*RTHT*RTHT*U+R* (R*UTTH+VTTH ) )/Sl  0296 

ELC= (RTTH*S1 +2.*RTHT*S2-R*S1 J/FMAG  0297 

FLC  =0.  0298 

GLC  =  — R* (SI *ST2  +  S2*ST*CT  J/FMAG  0299 

RETURN  0300 

END  0301 

SIBFTr  FSPDT.  LIST  0302 

SUBROUTINE  FSPDT ( RSP  »  THSP  «  FNT . THTSi THTFiOTHT . PH  I  « V X t  VY . VZ  )  0303 

C  INPUT  TH  TS » THTF iDTHT iPHI  tVXtVYtVZ  0304 

COMMON  AR1 « AR3.BR1 « AR9.AR1 Ot AR1 1 .WAVE  0305 

C  SEARCH  FOR  LARGEST  SCALAR  PRODUCT  IN  SECTION  0306 

FND=ABS( (THTF-THTS ) /DTHT )  0307 

ND=FND+ 1 •  0308 

FNT  =  0 •  0309 

RSP=0.  0310 

THSP=0 .  0311 

DO  1 00  I = 1 , ND. 1  0312 

F I  =  I  —  1  0313 

THT=F 1 *DTHT+THTS  0314 

R  =  R AD ( THT  )  0315 

CALL  FNORM(FNX.FNY,FNZ, R.THT. PHI )  0316 

FNE=VX*FNX+VY*FNY+VZ*FNZ  0317 

FNE=ABS (FNE )  0318 

1F(FNT.GT.FNF )  GO  TO  11  0319 
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FNTsABS(FNE)  0320 

RSP  =  R  0321 

THSP=THT  0322 

11  CONTINUE  0323 

100  CONTINUE  0324 

RETURN  0325 

END  0326 

SIBFTC  FINT.  LIST  0327 

SUBROUTINE  FINTISSS.FCTI . FLL  .  FUL . ERRR . NX )  0328 

INPUT  FCTI .FLL. FUL. ERRR. NX  0329 

EXTERNAL  DECLARATION  FOR  FUNCTION  FCTI  REQUIRED  0330 

COMPLEX  SSS. SS.FSS. FCTI .TRAP. TRAZ. SIMP, SIMZ.FNCP.FNCZ  0331 

FN=  NX  0332 

DEL=  <  FUL-FLL ) /FN  0333 

SSS= ( 0 . , 0  .  )  0334 

ERR=0 . 0 1 *ERRR/FN  0335 

A=FLL  0336 

DO  40  NNX= 1 . NX , 1  0337 

MXX  =  0  0338 

B=A+DEL  0339 

SS= ( 0 . , 0. )  0340 

MX=2  0341 

DX=DEL/2 •  0342 

LX= 1  0343 

X  =  A  0344 

GO  TO  15  0345 

5  TRAZ=DX*SS  0346 

MX=1  0347 

LX= 1  0348 

DX=DEL  0349 

10  SS=0.  0350 

LX=LX+1  0351 

DX=0.5»DX  0352 

X=A+DX  0353 

15  DO  20  I  X=  1  * MX » 1  0354 

FSS=FCT I(X)  0355 

SS=  SS  +  FSS  0356 

20  X  =X  +  2 • *DX  0357 

IF(LX.EQ.I)  GO  TO  5  0358 

MX=2*MX  0359 

TRAP=0 .5»TRAZ+DX#SS  0360 

DIF=CABS(TRAP-TRAZ)  0361 

I F (D IF.GE.D IP )  MXX  =  MXX+ 1  0362 

D I P=D I F  0363 

S I MP  = ( 4.*TRAP— TRAZ )/3.  0364 

FNCP= < 16.*SIMP-SIMZ)/15.  0365 

ER=CABS ( 1 .-FNCZ/FNCP )  0366 

TRAZ=TRAP  0367 

S I MZ  =S I MP  0368 

FNCZ  =FNCP  0369 

IFILX.LT. 4)  GO  TO  10  0370 

I F ( MXX . GT . 4  )  GO  TO  30  0371 

IF (ER.GT .ERR )  GO  TO  10  0372 

30  SSS=SSS+FNCP  0373 

40  A  =  A  +  DEL  0374 

50  CONTINUE  0375 

RETURN  0376 

END  0377 

SIBFTC  FFNRM.  LIST  0378 

SUBROUTINE  FNORM ( FNVX , FNV Y * FNVZ » R » THT » PH  I  )  0379 

C  INPUT  R, THT. PHI  0380 

COMMON  ARl , AR3.BR I , AR9, AR1 0. AR1 1 .WAVE  0381 

ST=SIN(THT)  0382 

CT  =  COS ( THT )  0383 
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SP=  S I N ( PH I  )  0304 

CP  =  COS ( PH  I  )  0385 

U=ARl*ST*ST+AR3*CT*CT+BRl*ST*CT  0386 

V=AR9*CT+AR10*ST  0307 

UTH=2 .* ( AR1-AR3 )*ST*CT+BRl * (CT*CT-ST*ST  )  0308 

VTH=-AR9*ST+AR1 0*CT  0389 

F 1 =2 • *R*U+V  0390 

F2=R*UTH+VTH  0391 

FX=ST*CP*F 1 +CT*CP*F2  0392 

F Y  =  ST *SP*F 1 +CT *SP*F  2  0393 

FZ  =  CT  *F 1 —ST  *F  2  0394 

FN=SQRT (FX*FX+FY*FV+FZ*FZ )  0395 

FNVX=FX/FN  0396 

FNVY=FY/FN  0397 

FNVZ  =  FZ /FN  0390 

RETURN  0399 

END  0400 

S1BFTC  RADD  LIST  0401 

FUNCTION  RAD(THT)  '  0402 

C  CALCULATE  R  0403 

COMMON  ARl . AR3. BR 1 » AR9, ARl 0 t ARl 1 . WAVE  0404 

ST  =  S 1 N ( THT )  0405 

CT=COS(THT)  0406 

C 1 =  AR1 *ST*ST  +  AR3*CT*CT+BR1  *ST*CT  0407 

C2=AR9*CT+AR 1 0*ST  0408 

C3= AR 1 1  0409 

IF (Cl .EQ.O. )GO  TO  11  0410 

ARG=SQRT (C2*C2-4.*C1*C3)  0411 

R =  (  — C2 +ARG ) / ( 2 • *C 1 )  0412 

R2= ( -C2-ARG ) / (2«*C 1  )  04  13 

IF (R2.LT.R.AND.R2.GT.0. )  R=R2  0414 

GO  TO  12  0415 

11  R  =  — C3/C2  0416 

12  R AD=R  0417 

RETURN  0418 

END  0419 

SIBFTC  FCOM.  LIST  0420 

SUBROUTINE  FCOMM(I)  0421 

COMMONRA 1 «RA3«RB 1 «RA9,RA 1 0 ,RA 1 1 « WAVE/D AT A/AR 1(20) « AR3 (20 ) iBR 1  (  20 ) «  0422 

CAR9(20),AR10(20) »AR1 1 (20)  0423 

RA 1 =AR 1 ( 1 )  0424 

R A3= AR3 ( I )  0425 

RB 1 =BR 1 ( 1 )  0426 

R A9= AR9 ( 1 )  0427 

RA 1 0= AR 10(1)  0428 

R A  1  1 =ARl  1(1)  0429 

RETURN  0430 

END  0431 

$DAT  A  0432 
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APPENDIX  II 

THE  COMPUTER  FLOW  DIAGRAM 


The  flow  diagram  presented  in  Ref.  1  is  presented  here. 
This  diagram  shows  the  sequence  in  which  the  operations  are 
performed. 
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